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ABSTRACT 


The nested mesoscale modelling system for air 
pollutants, developed under Phase I of this project, has been 
changed, to cover an area of approximately 750 x 600 km. The 
horizontal resolution is 20 km. The new domain includes Southern 


Ontario and the neighbouring portions of adjacent states. 


The model has been applied to simulate a high ozone 
episode in Southern Ontario. Based on the test runs, the 
convective boundary layer parameterization in Gesima was improved 
and optimal values of some input parameters were selected. 
Analyses were made on the areal distribution of ozone, NO,, and 
vocs and the time series of these species in selected locations: 
Guelph, Toronto and Oshawa. The results are in fairly good 
agreement with observations, which is a marked improvement over 
the results obtained in Phase I with respect to the predicted 


peak ozone levels. 
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EXECUTIVE SUMMARY 


Two existing numerical models: a German mesoscale 
meteorological model, Gesima and a long range transport model, 
Acid Deposition and Oxidant Model (ADOM), have been enhanced, 
adapted and interfaced under Phase I of this project to forma 
nested mesoscale modelling system for simulations of the 
transport, chemical transformations and deposition of atmospheric 
pollutants. Based on the tests of the model undertaken under 
Phase II, and described in this report, several significant 
improvements of the system were made. The domain of simulation 
has been extended from the domain used in Phase I to cover most 
of Southern Ontario as well as adjacent areas in the U.S. The 
horizontal resolution of 20 km used in the model, coupled with 
the appropriate topography, land use and local emissions data 
would allow finer local features of ambient pollution and 


deposition to be simulated. 


Two series of tests of the system were undertaken. The 
test runs of ADOM, showed that the ozone concentration is very 
sensitive to the inflow boundary concentrations. In order to 
decrease this sensitivity to boundary inputs, the domain has been 
extended from that used in Phase I. The second series of test 
runs focused on selecting the best set of input parameters for 
Gesima. The results of these tests also led to the replacement of 
the convective boundary layer parameterization in Gesima and some 


other modifications of the model. 


In order validate the system, a Simulation of a high 
ozone episode in Southern Ontario has been performed. The test 
run was carried out for a period of 8 days, starting on July 30 
1988. The analysis of the results of this run focuses on the 
areal distribution of ozone, as well as nitric oxides and 
hydrocarbons, and their diurnal variations in selected locations. 
The horizontal cross-sections of the ozone, NO, and VOC 
concentrations, presented in Figs. 14 to 28 show pronounced 
mesoscale variability, which cannot be simulated using regional 
scale air pollution models such as ADOM with horizontal 
resolution of 127 km as it is used to simulate over eastern North 
America now. The predicted levels of ozone at three Ontario 
locations (Guelph, Toronto and Oshawa) were compared to the 
available observational data. The time series (see Figs. 31 - 34) 
show a fairly good agreement with observations, especially near 
‘ground level. There is a marked improvement over the results 
obtained in Phase I, with respect to the predicted peak ozone 


levels. 
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1. INTRODUCTION 


This report describes the second phase of the development 
and testing of a nested grid, mesoscale air pollution modelling 
system, carried out under a contract with the Ontario Ministry of 
the Environment from July 1991 to June 1992. The test runs of 
the system are made for a high ozone episode over Southern 
Ontario in August 1988. 


The system is capable of simulating the emission, transport, 
chemical transformation and deposition of a wide range of air 
pollutants. Horizontal resolutions of 5 to 20 km can be used. 

The size of the domain of simulation can vary, depending on 
resolution and the available computer resources, from about 

100 x 100 km to 700 x 700 km or more. The test runs of the 
system for the Toronto area were carried out with horizontal 
resolution of 20 km in a domain of 340 x 340 km in Phase I. That 
domain has been enlarged in this phase to 750 x 600 km without 
changing the resolution. Thus, the present domain of the system 
covers most of Southern Ontario. 


The resolution and the size of the domain place this system 
in the mesoscale, i.e. between the urban scale air quality models 
and the large scale models of long range transport of atmospheric 
pollutants (LRTAP). The mesoscale modelling capability is 
necessary to address the urban NO,/VOC issue with a view to 
understanding and eventually to controlling the ozone precursor 
emissions and to quantify the subgrid scale variability of 
pollutant concentrations simulated by the LRTAP models. 


The approach was to adapt, enhance and interface two models: 
the large scale Acid Deposition and Oxidant Model (ADOM), and a 
mesoscale meteorological model ’Gesima’. Gesima and the modified 
mesoscale version of ADOM are then ‘nested’ in the large scale 
ADOM domain. 


The modelling system developed in this study consists of 
three main parts: the large scale ADOM, Gesima, and the 
mesoscale version of ADOM. The programs interfacing these three 
parts are also an important part of the system. The interactions 
between the system components are schematically represented in 
Figure 1. The initial and boundary conditions for the mesoscale 
ADOM are derived from the large scale ADOM output files made 
available by the Ontario Ministry of the Environment. 


The original large scale version of ADOM uses meteorological 
fields generated by the Canadian Meteorological Centre (CMC) 
spectral. diagnostic model. These data, having horizontal 
resolution of 381 km, are interpolated and gridded to the 127 km 
resolution domain by the ADOM meteorological preprocessor. The 


preprocessor also enhances the vertical resolution of the CMC 
fields by running a one-dimensional boundary layer model between 
the ground and the 2 km level at each grid point. For details see 
Scholtzeetwall 986). 


The mesoscale ADOM requires three-dimensional fields of 
meteorological parameters like horizontal and vertical wind 
velocity, temperature, humidity, turbulent diffusion parameters 
etc., with horizontal resolution of 20 km. Such fields are, in 
this project, generated by Gesima. The role of the large scale 
ADOM is to provide the initial and boundary conditions for the 
mesoscale simulations. Most of the development work of the system 
was completed under Phase I of this project. Details are 
described in the previous report (Niewiadomski and Shenfeld, 
1991). The approach of the project, major components of the 
system and the results of Phase I were also summarized in a paper 
presented at the Technology Transfer Conference in November 1991 
(Niewiadomski, 1991, included in this report as Appendix A) and 
in a poster presented at the Joint International Conference on 
Atmospheric Chemistry organized by the Canadian Institute for 
Research in Atmospheric Chemistry and the Ontario Section of AWMA 
in January 1992. 


Phase II of this project focused on testing, fine tuning and 
improving the system. Two types of tests were made, studying the 
response of ADOM to changes in the boundary conditions and point 
source emissions within the domain and the response of Gesima to 
changing the "nudging" coefficient. The optimal choice of other 
input parameters of the model was also determined. These test 
runs and their results are described in Sections 2 and 3. The 
tests led to a modification of the boundary layer module of 
Gesima which is described in Section 4. Section 5 presents some 
results of the first "production runs" of the system performed 
for the Ministry of the Environment in May and June 1992 as a 
part of the emission trading scenarios project. The adopted set 
of input parameters is provided in Appendix B. Appendix C 
describes in detail the modifications of the Gesima subroutines 
undertaken during Phase II. Short descriptions of all Gesima 
subroutines are included in Appendix D. 


A significant amount of time was also devoted to the 
technical problems resulting from changes in hardware and 
software due to the transfer of the system in April 1992 from the 
Cray computer of the Ontario Centre for Large Scale Computation 
to the Cray and SX-3 computers at the Environment Canada Centre 
Informatique de Dorval. The system is now fully operational at 
C.I.D. Its Gesima part operates on the new SX-3 supercomputer, 
while ADOM, containing some Cray-specific library subroutines, is 
currently run on the C.I.D Cray. 


2. TORONTO AREA TESTS 


At the beginning of this project several test runs of the 
system were performed using the same domain of simulations as in 
Phase I. These tests focused on the relative importance of 
boundary conditions and local point source emissions on ozone 
concentrations simulated by the mesoscale version of ADOM. 


The tests period covered August 1 to August 4, 1988. The 
mesoscale ADOM was run using the meteorological fields generated 
by Gesima during Phase I. For details on Gesima runs and the map 
of the domain, see the Phase I report (Niewiadomski and Shenfeld, 
1991). The output files of large scale ADOM, which provided the 
initial and boundary conditions for these tests, were produced by 
an improved version of ADOM that included a mass conservation 
module. 


The mesoscale ADOM was run with artificially increased and 
decreased boundary concentrations of ozone and NO,. Tests were 
carried out including and excluding two large NO, sources in the 
area: Nanticoke and Lakeview generating stations. Selected 
results of these tests are presented in Figures 2 - 7. 


The comparison of Figures 2 and 3 shows that a 100% increase 
in boundary ozone concentrations resulted at 8:00 p.m. EDT in an 
increase of ozone by 50 to 80% throughout the domain, without 
radically changing the areal distribution pattern. Additional 
increase of the NO, boundary values (see Figure 4) does not 
introduce significant changes. For the Toronto location (see 
Figure 6), the afternoon O, maxima increased by about 80% and the 
effect of an additional increase of boundary NO, values was 
negligible. Reducing the boundary concentrations by 90% resulted 
in the reduction of maximum 0, by 60% to about 25 ppb. Thus, the 
emissions from inside the domain contributed only to 40% of the 
simulated base case maximum Toronto ozone concentration. The 
contributions of single sources is much lower, as can be seen 
from Figure 7 which shows O, time series for the case with 
Lakeview emissions excluded. 


The dependence on boundary conditions, which are derived 
from the large scale ADOM results and thus do not have sufficient 
resolution, brings into question the validity of using such a 
small domain for mesoscale simulations. It was decided that 
instead of running the model in smaller domains separately e.g. 
for the Toronto and Windsor areas, the domain should be large 
enough not only to cover most areas of interest but also to move 
the boundaries of the domain far from those areas. As can be 
seen from Figure 8, the new larger domain fulfils the latter 
condition much better for Toronto than for Windsor. 


ah 


3. SOUTHERN ONTARIO TESTS 


3.1 Domain of Simulations and the Set-Up of the Experiment 


The new large domain of simulation is shown in Figure 8. It 
covers 6 x 5 grid cells of the large scale ADOM domain (with x 
and y coordinates 13 to 18 and 15 to 19, respectively). To 
obtain the horizontal resolution of about 20 km this area has 
been divided into 37 x 29 cells. Note that due to the 
peculiarities of the graphic software used, the grid shown in 
Figure 8 and in the subsequent figures does not contain zero 
numbered lines, so the grid lines in the figures correspond to 
the centres of the grid cells rather then to their borders. 


As in Phase I, the basic vertical structure of the large 
scale ADOM, with 12 variable depth layers, has been retained for 
the mesoscale version of the model. The tops of the ADOM layers 
aremlocated at 5652, 1135-8/0250.7, 416-3, 655-3, 100050, 1497 22, 
2214.5, 3249.2, 4741.6, 6894.5, and 10000.0 m above the ground 
level. As in Phase I, Gesima uses 16 vertical layers, with each 
of the 4 top layers of ADOM corresponding to two Gesima’s layers. 
Thus the tops of Gesima layers lie at 56.2, 135.8, 250.7, 416.3, 
655.3, 1000.0, 1497.2, 2214.5, 2731.8, 3249.2, 3995.4, 4741.6, 
5818.1, 6894.5, 8447.3, and 10000.0 m. The Gesima/ADOM interface 
program averages the data from top layers of Gesima and assigns 
them to the appropriate layers of ADOM. Note that the heights of 
the layer tops given above are relative to the average 
topographic height of the grid cell. 


Similar to Phase I, the well documented high ozone episode 
in August 1988 has been selected for simulations, but the period 
has been extended from 4 to 8 days, starting on July 30 1988. 

The large scale ADOM output files were obtained by a new 
simulation on a reduced (33 x 21) horizontal domain. Twelve 
southern rows of the original grid were removed and higher 
concentrations ozone were imposed on the new southern boundary, 
where ozone levels had been underpredicted by the ADOM model in 
previous runs. These higher values were derived from hourly 
surface observations interpolated throughout the boundary layer 
to match predictions of ADOM run on the full domain (33 x 33 grid 
cells. The boundary conditions for other species were derived 
from a former simulation on a full grid. The southern boundary of 
the mesoscale domain is located 254 km north of the south 
boundary of the 33 x 21 grid domain of the large scale ADOM. 


Some input data files of Gesima and mesoscale ADOM have been 
prepared in a different way than in Phase I. The Ministry of the 
Environment made available the topography and land use data files 
for the mesoscale domain, so these files were used directly as 
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input to Gesima and to prepare the geophysical input file of the 
mesoscale ADOM. This made redundant some interpolation 
algorithms in Gesima. Furthermore, since the land use classes 
used in Phase I in Gesima differed from those of ADOM, the Gesima 
classes have been redefined so that the same input data sets can 
be used for both models. For details see Section C.2 of the 
Appendix C. 


3.2 Results 


The tests described in this section focus on the sensitivity 
of the Gesima results to the value of the nudging coefficient. 
For the discussion of the nudging terms see Section 2.4 of the 
Phase I report. Analysis of the results of these runs led not 
only to the selection of proper nudging coefficients, but also to 
several improvements of the system and to better selection of 
other input data. The most important improvement was the 
modification of the boundary layer module of Gesima, described in 
detail in Section 4. Others included minor modifications of the 
surface layer module, and changes to the radiation scheme (for 
details see Appendix C). The values of the input parameters 
selected as a result of these tests are listed in Appendix B. 


Tests also showed that the boundary conditions scheme for 
the cloud variables in Gesima (see Section A1.5 of the Phase I 
report) introduced some errors to the results. The horizontal 
resolution of 20 km was also too coarse to properly simulate the 
convective clouds. Since the clouds do not play an important role 
in the episode studied, it was decided to exclude the whole cloud 
module, until its boundary condition scheme can be fully tested 
and improved. The effects of clouds on incoming solar and 
atmospheric radiation in the surface layer are parameterized 
according to the interpolated large scale cloudiness data. For 
details see Appendix C. 


Several test runs with different nudging coefficients were 
performed. As expected, increasing the value of the nudging 
coefficient led to greater suppression of mesoscale effects, and 
consequently greater uniformity over the domain while decreasing 
the nudging coefficient resulted in the mesoscale fields becoming 
inconsistent with the large scale meteorological situation. 

Based on these tests, a vertically uniform nudging coefficient of 
0.0001 per second was chosen for all variables. The surface wind 
fields simulated with this coefficient for August 2, 1988 are 
compared to the winds i.e. interpolated from the large scale data 
in Figures 9 - 13. As can be seen the mesoscale winds, while 
retaining some local features follow the main large scale flow. 


Even so, it takes several hours for any significant change in the 
large scale situation to be reflected in the simulated mesoscale 
fields. 


4. NEW CONVECTIVE BOUNDARY LAYER SCHEME FOR GESIMA 


The original Gesima determined the turbulent diffusivity 
coefficients for momentum and heat at the height z according to 
the following formulae: 


Ke, j= (mt /02) IS eae 00707210 (1) 
Kym = m? ($?2-15B) * if 00/dz< 0 (2) 
Ko Kun D (3) 


where ®. is a stability function computed according to Bussinger 
et al. (1971), m/=*k2 7 G@: +5 kz/-X.)) -ist@the maxing length 
computed according to Blackadar (1962) and Mellor and Yamada 
(1974). S = ((9u/@z)* + (3v/8z)‘)" and B = (gf) (49/9z) denote 
the wind shear and buoyancy parameter, respectively, 

while k = 0.4 is the von Karman’s constant. D is a coefficient 
based on the Richardson’s number and à. is a function of the 
Coriolis parameter and the mean geostrophic wind. 


The model tests showed that for resolution of 20 km the 
above formulation, which was originally designed to be run with 
horizontal grid sizes of 2 to 5 km, led to underestimated 
diffusivities in cases of convective conditions. The mixed layer 
extended up only several hundred meters or less, even on a sunny 
summer afternoon when the mixing height could be up to 2000 


meters. 


Equation (2) for K. during convective conditions was changed 
based on the approach of Hass (1991) to the following: 


Kae = 2k swe Zz (2S = 2/25) (4) 
where Z, is the mixed layer height. 
The convective velocity scale (w.) is computed from the formula: 
we = (GHZ /(T¢0))? (5) 


where and c, are density and specific heat of air and T denotes 
the avérage temperature in the mixed layer. The surface heat flux 


7) 


(H) is computed in the surface layer module of Gesima. Since 
the mixed layer height (Z) is not directly available in Gesima, 
its value is interpolated from the low resolution input data of 
the large scale ADOM. 


The criterion for applying the convective formulation is now 
based on the sign of H so that equation (4) is used for z < 2, 
only if the surface heat flux is directed upwards. In all other 
cases equation (1) applies as before. The diffusivity 
coefficient for heat is computed according to equation (3) as in 
the original Gesima. 


The difference between the new and original formulations 
is illustrated below in Table 1, showing vertical profiles of 
turbulent diffusion coefficient at 4.00 p.m. EDT on August 1, 
1988 for a grid cell NE of Toronto. Z, for that case was equal to 
1130 m, and w. to 0.50 m/s. 


Table 1. Sample diffusivity profiles 


original formulation |modified formulation 
ro 














136 


56 0.328E+02 
251 
416 


655 


1000 
1497 
2214 
2731 


rR 


3995 
4742 
5818 
6894 
8447 
10000 


= 
N 


FH 


br IR 
un 
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3249 0.366E-02 0.369E-02 | 0.366E-02 | 0.369E-02 


= 
aS 
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se EMISSION TRADING RUNS 
5.1 Introduction 


After the completion of the Gesima tests described in 
Section 3, several runs of the system were performed as a part of 
the emission trading project of the Ontario Ministry of the 
Environment. The meteorological fields were produced by Gesima 
for the period July 30 to August 6, 1988 in the domain described 
in Section 3.1 and shown in Figure 8. Mesoscale ADOM was run for 
the same period, using emission data provided by O.M.E., 
representing 4 different emission scenarios. 


The results for various scenarios will be described in a 
separate report. The ozone distribution for the base case, which 
has the input emissions most closely corresponding to the actual 
emissions in 1988, is discussed below. 


5.2 Areal Distribution of pollutants 


Contours of the level 1 simulated ozone concentration at 
4:00 p.m. EDT (i.e. close to the daily ozone peaks) on each day 
of the period are presented in Figs. 14 to 21. Figure 22 shows 
the ozone concentration at 4 a.m. on August 2 1988 and Figs. 23 - 
28 the concentrations of NO, NO, and VOCs at 4 a.m. and 4 p.m. of 
that day. The VOC levels in Figs. 27 - 28 (and Figs. 29 - 30 as 
well) represent combined concentrations of 11 ADOM species, 
computed according to a formula provided by O.M.E. 


As expected, the areas of highest ozone concentration are 
located downwind of major urban areas, Detroit/Windsor, Cleveland 
and Toronto. At night these areas have the lowest oO, 
concentrations due to scavenging by NO,. The higher ozone 
concentrations over the Great Lakes are caused by the lower 
diffusivity and dry deposition over water than over land. 


5.3 Diurnal Variations and Comparison with Observations 


The time evolution of concentrations of ozone, NO, NO, and 
vVoCs modelled at level 1 of the model on August 2, 1988 are shown 
for the Toronto and Guelph cells in Figs. 29 and 30, 
respectively. At both locations NO, and VOCs peak at the same 
time: 9 a.m. EDT. The afternoon peak of ozone occurs in Toronto 
at 4:00 p.m, one hour later than in Guelph. The concentrations of 
NO and VOCs are an order of magnitude greater in Toronto than in 
Guelph. Toronto has twice the Guelph NO, concentration. The ozone 
peak in Toronto is, however, only slightly higher than in Guelph, 


indicating the importance of the ratio of VOC to NO, rather than 
their individual magnitudes in the development of ozone. 


The plots of the time evolution of the modelled ozone 
concentration in selected locations are given in Figs. 31 to 34. 
The results for the model layers 1 and 4 of the Toronto cell 
(Figs. 31, 32) are given for the whole period of simulation, 
while those for the lowest level of the grid cells containing 
Guelph (Fig. 33) and Oshawa (Fig. 34) cover only the middle part 
of the period (August 1 to 4). Note that similar plots were 
presented for the same locations in the report of Phase I 
(Niewiadomski and Shenfeld, 1991). 


In order to validate the ozone levels as computed by the 
ADOM/Gesima approach, the simulated concentrations are compared 
to the levels measured from August 1 to 4, 1988, in one hour 
intervals, at the Air Quality Index (AQI) stations at the heart 
of downtown Toronto (Fig. 31) and at the top of CN Tower (Fig. 
32). Similar measurements from the AQI stations in Guelph and 
Oshawa are included in Figs. 33 and 34. 


It should be noted that while the measurements were taken at 
specified points, the modelled concentrations represent averages 
from a mesoscale ADOM cell of the size 20 x 20 km and a depth of 
56 m for level 1 and a depth of 165 m for level 4. Even so, the 
modelled and observed ozone concentrations are in generally good 
agreement. Although the model does not predict the highest 
measured concentrations (e.g. 112 ppb observed in Toronto on 
August 2, 1988 at 14:00 EDT versus 75 ppb modelled), such 
discrepancies result, at least partly, from the incompatibility 
between point measurements and the model’s volume averaged 
concentrations. 


The times of maximum and minimum modelled concentrations 
coincide with those observed at the surface in all three 
locations (except for the secondary peaks in observed 
concentrations at the early morning of August 4, 1988, in Toronto 
and Oshawa, which are not reflected in the model results and a 
secondary peak in model results at the early morning of August 2 
in Guelph which is not reflected in measurements). In general, 
the modelled minima tend to have lower concentrations and last 
longer than the observed ones. 


The upper level data, both modelled and observed, have much 
less regular diurnal variability (see Fig. 32). The CN Tower is 
close to the tall buildings in downtown Toronto, so the 
turbulence structure is rather complicated. This may explain the 
large fluctuation in the observed ozone concentrations. At the 
higher altitude there is no removal of ozone due to deposition 
and less scavenging due to the automobile emissions of nitric 
oxide than at ground level. The modelled concentrations generally 
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agree with the observed levels. The largest discrepancy occurs at 
6:00 EDT on August 4, 1988. Neither the peak in the model 
results at that time, nor the minimum in observed levels can be 
explained by meteorological or other input parameters. 


A remarkable improvement, especially for the lowest layer, 
can be seen when comparing Figs. 31 - 34 to corresponding results 
of Phase I (Figs. 113 - 116 of Niewiadomski and Shenfeld, 1991 
and Fig. 2 of Appendix A of this report). The increase in the 
modelled ozone concentrations during the afternoon peaks (which 
however still do not, in several cases, reach the observed 
values) can be attributed to improvements in the input emission 
files, southern boundary conditions of the large scale ADOM and 
the improvements of the diffusivity parameterization in Gesima. 
Improvement in the timing of the ozone peaks, especially in 
Guelph, results probably from increasing the model domain (which 
moved the western boundary of the domain much further from that 
station) and the mass conservation correction in the large scale 
ADOM, introduced after the completion of Phase I. 


6. SUMMARY AND CONCLUSIONS 


A mesoscale modelling system for simulating the transport, 
chemical transformations and deposition of atmospheric 
pollutants, designed, developed and preliminarily tested under 
Phase I of this project, has been extended to cover most of 
Southern Ontario. Knowledge gained from the test runs resulted in 
significant improvements of all components of the systen, 
especially the turbulent diffusion parameterization in Gesima, 
and better selection of adjustable input parameters. 


In May and June 1992 the first "production runs" of the 
system have been performed, as a part of the Ontario Ministry of 
the Environment emission trading project. 


Preliminary analyses of the simulation of the August 1988 
high ozone episode over Southern Ontario showed reasonably good 
agreement of the model results with observations and a marked 
improvement over the Phase I results. 


The system is capable of simulating the distribution of 
various air pollutants with the spatial resolution of 20 km or 
less. It can detect local effect of this scale, which is 
impossible to simulate with large scale models like the original 


version of ADOM. 


More comprehensive analysis of the results of several case 
studies, and further improvements of the system, including 
decreasing the grid size to 5 x 5 km, are planned for Phase III 
of this project to be carried out from July 1992 to June 1993. 
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Figure l. Components and scheme of operation of the nested mesoscale 
modelling system. 
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Figure 2: Simulated ozone concentrations (in ug per cubic meter) 
at 8:00 p.m. EDT on August 2, 1988, base case. 
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Figure 3: Simulated ozone concentrations (in ug per cubic meter) 
at 8:00 p.m. EDT on August 2, 1988, boundary concentrations 
of ozone increased by 100%. 
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Figure 4: Simulated ozone concentrations (in ug per cubic meter) 
at 8:00 p.m. EDT on August 2, 1988, boundary concentrations 
of ozone and NOx increased by 100%. 
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Figure 5: Simulated ozone concentrations (in ug per cubic meter) 
at 8:00 p.m. EDT on August 2, 1988, boundary concentrations 
of ozone and NOx increased by 100%, Lakeview and Nanticoke 
sources removed. 
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Figure 7. O base case + Lakeview, Nanticoke removed. 
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Figure 9: Surface wind fields at 00:00 GMT on August 2, 1988; 
top: interpolated from large scale data, 
botton: simulated by Gesima 
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Figure 12: Surface wind fields at 18:00 GMT on August PRISES; 
top: interpolated from large scale data, 
bottom: simulated by Gesima 
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APPENDIX A 


THE DEVELOPMENT OF A LONG RANGE TRANSPORT MODEL 
WITH A NESTED FINE RESOLUTION GRID 


by 
M. Niewiadomski, The MEP Company, Markham, Ontario L3R 9T2 


l. INTRODUCTION 


Two existing numerical models, GESIMA, a mesoscale, non- 
hydrostatic meteorological model developed at the GKSS 
Forschungszentrum in Geesthacht, Germany and the long range 
transport Acid Deposition and Oxidant Model (ADOM) have been 
enhanced, adapted and interfaced to obtain a nested modelling 
system for the simulation of the transport, chemical 
transformation and deposition of a wide range of atmospheric 
pollutants. 


The system described here is best suited to simulate the 
distribution of pollutants with a horizontal resolution of 5 to 
20 km. A test run of the system for the Toronto area was carried 
out in a domain of 340 x 340 km with horizontal resolution of 20 
km. A new series of simulations, with the same resolution, is 
currently being performed in a domain of 700 x 500 km covering 
most of Southern Ontario. 


The resolution and the size of the domain place this system 
in the mesoscale, i.e. between the urban scale air quality models 
and the large scale models of long range transport of atmospheric 
pollutants (LRTAP). A mesoscale modelling capability is 
necessary to address the urban NOx/VOC issue if an understanding 
and eventual controlling of the precursor emissions is to be 
achieved. This approach is also required and to quantify the 
subgrid scale variability of pollutant concentrations simulated 
by the LRTAP models. 


The modelling system developed in this study involves the 
interrelation of three main parts: the large scale ADOM, GESIMA, 
and the mesoscale version of ADOM. The programs interfacing 
these three parts are also an important part of the system. The 
interactions between these system components are schematically 
represented in Figure 1. 


The original large scale version of ADOM uses meteorological 
fields generated by the Canadian Meteorological Centre (CMC) 
spectral diagnostic model. The mesoscale ADOM requires 
three-dimensional fields of meteorological parameters like 
horizontal and vertical wind velocity, temperature, humidity, 
turbulent diffusion parameters etc., with much finer horizontal 
resolution. Such fields are provided in this project by GESIMA. 
The role of the large scale ADOM, is to provide the initial and 
boundary conditions for the mesoscale simulations. 
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For a brief description of the main features of the system 
components, see Niewiadomski (1990). More details can be found 
in the original papers on ADOM (Scire et al., 1986; Venkatram et 
al., 1988) and GESIMA (Kapitza, 1987; Kapitza and Eppel, 1987; 
Jacob et al., 1990 and Mengelkamp, 1991). 


2. THE NESTING SCHEME 


For simulations periods longer than, say, one day, a 
mesoscale model must receive information of the development of 
the modelled fields outside its domain. This is achieved by 
"nesting" the mesoscale model inside an outer model operating in 
a much larger domain but with a more coarse resolution. In the 
system described here the outer model for GESIMA is the CMC 
diagnostic model, which also provides the meteorological fields 
for the large scale ADOM. The mesoscale ADOM is nested in the 
large scale version operating on the same grid as the CMC model. 


In this system the mesoscale domain covers a specified 
number of cells of the large scale grid. The large scale data 
from these cells, and those adjacent to the mesoscale domain, can 
be interpolated in space and, if necessary, in time and then 
passed to the mesoscale model. 


The mesoscale ADOM receives information from the outer model 
through time dependent boundary conditions. Similar time 
dependent boundary conditions are used for some parameters in 
GESIMA, but for the main fields (wind, temperature and humidity) 
a technique of the Newtonian relaxation or "nudging" is used. 


The concept of nudging (Hoke and Anthes, 1976; Stauffer and 
Seaman, 1990; Seaman and Cole, 1991) involves the introduction of 
additional, artificial tendency terms to the prognostic 
equations, which force the solution of the mesoscale model toward 
the observed, or modelled, large scale conditions. This approach 
is used to combine the fine scale variability with the large 
scale tendencies which cannot be simulated by the mesoscale model 
itself, and seems to be more efficient and reasonable than 
introducing those tendencies by means of the time dependent 
boundary conditions alone. 


These artificial tendencies - the nudging terms, applied at 
each time step, are proportional to the difference between the 
values of large scale field at given grid point (interpolated 
from the results of the outer model) and the actual solution of 
the mesoscale model. The coefficient of proportionality - the 
nudging coefficient, can be a function of height and has values 
of the order of 0.0001/sec. 
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Note that the described above nesting scheme represents so 
called one-way nesting, i.e. the results of the mesoscale 
simulations do not influence the large scale models. This allows 
all components of the system to be run separately and to simulate 
several emission/pollution scenarios using the same 
meteorological fields. 


3. SIMULATION OF THE HIGH OZONE EPISODE (August 1 - 4, 1988) 


This section gives the details of a test run of the system. 
The episode for this study, August 1 - 4, 1988 was chosen because 
of the high ozone concentrations observed at that time in 
Southern Ontario, and the availability of the results of the 
large scale ADOM simulation and the observational data. 


Some results of this simulation concerning the time 
evolution of the ozone concentrations in Toronto and its areal 
distribution in the whole domain are presented and discussed in 
Section 4. 


The large scale ADOM is routinely run on a 33 x 33 grid 
based on a polar stereographic projection with nominal horizontal 
resolution of 127 km (true at 60° N). That domain covers most of 
eastern North America. The same grid is used in the CMC spectal 
model. 


In this simulation, the ‘nested’ grid of GESIMA and 
mesoscale ADOM covers 9 cells of the large grid, including the 
Toronto cell and the 8 cells around it. Thus, the mesoscale 
domain extends over an area of about 120,000 km with 
Metropolitan Toronto area located about 25 km south of its 
centre. The domain is divided into 17 x 17 mesoscale cells 
resulting in a horizontal resolution of about 20 km. 


The basic vertical structure of the large scale ADOM, with 
12 variable depth layers has been retained for the mesoscale 
version of the model. The same structure was adopted for the 
GESIMA simulations. The tops of the ADOM layers are located at 
D632, 135.0, 250.7, 416.3, 655.077 L000 .0,.1497 .2 4 2214.5, 3249 2, 
4741.6, 6894.5, and 10000.0 metres above the ground level. The 
vertical resolution of the upper layers of ADOM is too coarse for 
GESIMA, and therefore the 4 top layers of ADOM were split into 
two layers. In effect, 16 rather than 12 vertical layers are used 
in GESIMA. The GESIMA/ADOM interface averages the data from those 
split layers before using them in appropriate layers of ADOM. 


A4 


GESIMA was run for a period of 96 hours starting on 
August 1, 1988 at 00:00 GMT (July 31, 20:00 EDT) and a time step 
of 45 sec was used throughout this period. 


The large scale ADOM was not run under this project. The 
output files, from a previous simulation for the same period were 
used. These files together with the output of GESIMA, were 
processed by the interfacing programs and provided the input data 
for the mesoscale ADOM simulations. 


Since the mesoscale ADOM requires GESIMA output data, the 
ADOM run started one hour later then that of GESIMA, at 01:00 
GMT, on August 1, 1988 and continued for 95 hours. 


4. RESULTS 


The time evolution of the ozone concentration in the Toronto 
cell of the mesoscale model is compared with observations in 
Figure 2. 


The times of maximum and minimum modelled concentrations 
coincide with those observed at the surface in Toronto (except 
for a secondary peak in observed concentrations during the early 
morning of Aug 4 which is not reflected in the model results). 


The values of the modelled ozone concentrations are also in 
fairly good agreement with the observed ones, although the model 
does not predict the highest measured concentrations (e.g. 112 
ppb observed in Toronto on Aug. 2 at 14:00 EDT versus 68 ppb 
modelled). These discrepancies can be partly attributed to a 
significant area and depth represented by an ADOM cell, compared 
to the point measurements. Comparisons made for several other 
locations with available surface measurements of ozone 
concentrations yielded similar results. 


Examples of the areal distribution of ozone and nitrogen 
oxides are shown in Figure 3. While the NOx distribution is 
dominated by strong area and point source emissions in the 
Toronto area and the ozone is distributed more evenly, both 
fields display a pronounced horizontal variability impossible to 
simulate with a large scale - coarse resolution models. 


5. CONCLUSIONS 


A mesoscale modelling system for transport, chemical 
transformation and deposition of atmospheric pollutants has been 
designed, programmed and tested in one case study over Southern 
Ontario. 


The system is capable of simulating the distribution of 
various air pollutants with a spatial resolution of 20 km or 
less. It can detect local effects in this scale that are 
impossible to simulate with large scale models like the original 
version of ADOM. 


A preliminary analysis of the simulation of the early August 
1988 high ozone episode over Southern Ontario, showed reasonably 
good agreement of the model results with observations. 


A case study of the same episode using a slightly modified 
model, over a large domain covering most of Southern Ontario is 
currently underway. 
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Description of the Gesima input file 
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APPENDIX B 
Description of the Gesima input file (fort.15) 


The general input file is organized in namelists. The file is 
read by subroutine READD at the beginning of each job. Namelists 
HOTCOL, SWITCH and AEROSO are also read from MAIN. Only those 
three namelists are read in the case of a restart (continuation) 
run, thus only the parameters belonging to those namelists can 
change values after the beginning of an initial job. 


The basic structure of the input file and most of the input 
parameters originate from the original, GKSS version of the 
model. Some of them are redundant for the present application, 
but have been retained in the input file for the sake of 
compatibility with original model and therefore require a value 
to be provided in the input file. Other parameters may require 
specific 

values for this study. Such values are provided below, where 
applicable. For other parameters either a recommended value (for 
the 37 x 29 grid with 20 km resolution and the August 1988 
simulation) is given or an exemplary value from a recent run. 


Le Namelist HOTCOL 


This namelist contains the main control parameters of the model. 
For the regular simulations, the model should be run for one or 
more full hours, beginning on IHRO and ending at IHR1 on the same 
date. Other simulation time parameters are of secondary 
importance. The length of the time step DELTAT should be chosen 
so that an integer number of steps is made in one hour. The 
description of input parameters is given below. 


Initial/restart run indicator: 
IFHOTS 0 for an initial run 
for a continuation run 


oil 
| 


Intervals (in seconds) of calls to SAFETY and writing a restart 
file: 
INCHOT = 3600 (exemplary value) 


Time step in seconds: 
DELTAT = 45. (recommended value, see comments on 
DELTAT above) 


Intervals (in seconds) of calls to PRINTT and writing output 
files: 
INCOUT = 21600 (exemplary value) 
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Intervals (in seconds) of calls to MONITO and writing auxiliary 
output: 
IFMONI = 20000 (exemplary value) 


Maximum total time of simulation: 

TIMAX = 10000000 (any large number) 

Setting TIMAX to lower values is useful for short test runs. The 
simulation will stop after reaching TIMAX, even if IHR1 is not 
reached. 


INRAD = 2 (any value - not used in this 
version) 
IFTIM = 0 (required value - TIMSTP not 


called, for IFTIM = 1 variable 
timestep is internally computed). 


IFSTAT = 0 (required value, IFSTAT = 1 means 
‘static’ simulation) 
Filter parameters: 


FILX = 0.025 (recommended value) 

FILY = 0.025 (recommended value) 

FILZ = 0.00 (recommended value) 

Smoothing parameters 

SMCO = 0.0 (any value - not used in this 
version) 

SMCOM = 0.0 (any value - not used in this 

version) 


Geostrophic wind parameter: 

KGEO = 10 (recommended value, geostrophic 
wind is assumed as equal to the 
nudging wind at appropriate level 
for all levels above KGEO and to 
the nudging wind at KGEO below 
that level) 


Initial and final hour of simulations (GMT): 
IHRO 0 (exemplary value) 
IHR1 6 (exemplary value) 


First hour (GMT) for which the nudging data should be read from 
units 91 and 41 rather than 81 and 21: 
IHR91 = 24 (exemplary value) 
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IFRELA = -1 (Value required for a regular run, 
other values can be used for some 
special, mostly obsolete cases, 
for details see subroutines 
RDBC and RDBCI) 


IDAY ST (any value less then 10, other 
values can be used for some 
special, mostly obsolete cases, 
for details see subroutines RDBC 
and RDBCI) 


Intervals (in seconds) of calls to WRTAPE and writing ADOM output 
file: 
INTAP = 3600 (recommended value) 


Nudging coefficients for humidity, temperature and wind, NZ+2 
values each: 


QNUD = 18*1.E-4 (recommended values) 

TNUD = 18*1.E-4 (recommended values) 

UNUD = 18*1.E-4 (recommended values) 

Model version: 

VRSION = 'GES/ENE’ (exemplary value, not important in 


this application) 


Start of simulation date: 


ISTART = 880806 (exemplary value) 
Julian time (day, GMT hour) of the start of simulation 
IJULS = 21900 (exemplary value) 
End of simulation date: 
IIEND = 880806 (exemplary value) 
ADOM output file append switch: 
ITAPPE = A0 if a new ADOM output file (fort.32) 
is to be written. 
= 1 If ADOM output is to be appended 


to the results of a previous run. 


2% Namelist GEOMET 


This namelist contains the grid definition and some auxiliary 
parameters. Only the vertical grid data are meaningful for this 
application. The vertical layers correspond to those of ADOM, 
with 4 uppermost layers split in two to enhance the resolution. 
The additional level at the top (NZ+1) is at the same distance 
from the top as the level NZ-1). 
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Horizontal grid (NX+3 values in x direction and NY+3 in y 
direction): 


XE = 40*0. (any values, grid redefined in 
subroutine READD) 
YN =,32%0.. (any values, grid redefined in 


subroutine READD) 


Vertical grid (tops of each layer in m, NZP+2 values), required 

values: 

ZW =100,56-1983,135:8103,12250:6563 
416-3205,655-2896/;241000., 1497.2413, 2214.5072, 
2731583146, 321945593995. 3912, 474156265, 
5818.0661, 6894.5056, 447.2528, 10000., 


11552..7472 
Height of the bottom and top of the domain (m): 
HBOT = 0. (required value) 
HTOP = 10000. (required value) 
RTOP = 25000. (any value - not used in this 
version) 


Height of topography (in m, (NX+3)*(NY+3) values): 
ZGNE = 1280*0. (any values, topography redefined 
in subroutine RTOPO) 


Land use category ae gna (NY+2) values: 
LWTYP 1209*8 (any values - not used in this 
version) 


3. Namelist PARMET 


This namelist contains mostly physical parameters, which do not 
change during the simulation. 


Average latitude of the domain (deg): 
PHI = 44, (exemplary value) 


Start of simulation in days from the beginning of spring 
DAY = 130. (exemplary value) 


Mean cloud coverage: 
COV = 0.0 (any value - not used in this 
version) 
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Sea surface temperature: 


TNSEA = 286. (any value - not used in this 
version) 

Pressure solver parameters: 

ITMAX = 200 (recommended value) 

RESMAX = 2.E-04 (recommended value) 

RESM2 = 2.E-2 (recommended value) 

RESFAC = 4. (recommended value) 

Geostrophic wind components (m/s): 

UGEO = 4.0 (any value - not used in this 
version) 

VGEO = -2.0 (any value - not used in this 
version) 

Reference potential temperature profile (NZ+2) values: 

THETO = 18*290. (recommended values) 

Reference sea level pressure (Pa): 

PSEAL = 1.E5 (recommended value) 

Reference sea level temperature (K): . 

TO = 290. (recommended value) 

Reference sea level potential temperature (K): 

THSEAL = 290. (recommended value) 

Reference temperature profile parameter: 

GAMMA = 0.6 (recommended value) 

Parameters for test Mer runs - see subroutine READD) 

AMPL = 200. (any value - not used in this 
version) 

AMPL2 = 300., (any value - not used in this 
version) 

WIDTH = 10000., (any value - not used in this 
version) 

WIDTH2 = 25000., (any value - not used in this 

version) 

XCENT = 75000., (any value - not used in this 
version) 

YCENT = 150000., (any value - not used in this 
version) 

Height of the damping layer 

ZDAMP = 6894.5056 (any value, damping layer defined 


by IFDAMP) 
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Damping coefficient: 

DAMFAC = 0.024 (recommended value) 

A small number: 

SMALL = 1.E-30 (recommended value) 

Implicit diffusion parameter: 

TDIF = 1.50 (recommended value) 

Tracer parameters, 3 values each: 

XTRAC = 1000., 5000, 10000. (any values - not used in this 
version) 

YTRAC = 1000., 3000., 3000. (any values - not used in this 
version) 

ZTRAC = WO 5 Ade, UHOwe (any values - not used in this 
version) 


Permutation parameter, 24 values (required values): 


4. Namelist SWITCH 


This namelist contains switch variables which activate or 
deactivate various options of the model. Usually only one choice 
is possible for this application. The purpose of each switch is 
given below along with the settings required for this study. For 
other possible values of the switches and their meanings see the 
appropriate subroutines. 


Note that even that this namelist is read also at the beginning 
of restart jobs, and the setting of switches can be changed, some 
of such changes may introduce inconsistencies in the model and 
lead to the execution time errors. 


Inflow and outflow boundary condition switches for temperature, 
wind, vertical velocity and humidity: 


IFTINF = 1 (required value) 
IFTOUT =a (required value) 
IFUINF =-1 (required value) 
IFUOUT =-1 (required value) 
IFWINF ele (required value) 
IFWOUT =1 (required value) 
IFQINF =1 (required value) 
IFQOUT = 1 (required value) 


Nudging switches for temperature, humidity and wind (1: ON, 


= 1 (required value) 
IFNUQ =) d (required value) 
= 1 (required value) 


4h Æ 
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Velocity initialization switches: 


IFINUW = i (required value) 
IFINFB = 1 (required value) 
IFINFL = 1 (required value) 


Thicknes of the _Gamping layer (IFDAMT uppermost levels): 
IFDAMT 4 (recommended value) 


Surface wetness switch: 


IFSWET = 0 (required value) 
Land use data switch: 

IFLWTY = 0 (required value) 
Temperature initialization switches: 

IFINTP = 0 (required value) 
IFTETA = 1 (recommended value) 


Reference temperature profile switch: 


IFPROF = 2 (required value) 
Filter switch (0,3 or 5 point filtering): 
IFFILT = 5 (required value) 
Topography Cyne. switch: 

IFTOPO 0 (required value) 
Turbulence parameterization switch: 

IFDIFF = 2 (required value) 
Hydrostatic pressure switches: 

IFHYUP =1 (required value) 
IFHYTO = 1 (required value) 
IFHYGR = 1 (required value) 


Auxiliary Qutpur switch: 
IFCHEC 0 (recommended value) 


Coriolis force switch: 
IFCORI = 0 (recommended value) 


Top and bottom boundary condition switches: 


IFTBOT eae) (required value) 
IFTTOP = 0 (required value) 
IFUBOT = 3 (required value) 
IFUTOP = 0 (required value) 
IFCBOT = 3 (required value) 
IFCTOP = 1 (required value) 


Humidity simulation switch: 
IFHUMI a (required value) 
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Cloud module switch: 
IFCLOU = 0 (required value) 
Surface module switch: 
IFSURF = 1 (required value) 
Passive tracer switch: 
IFPASS = 0 (required value) 
Longwave radiation switch: 
IFRAD = 0 (required value) 


Unit numbers for output files written by subroutine PRINTT, 

22 values, (recommended values): 

IFOUTP = D 227 53, 047 55), 90,7 57, 58, 59, 60, 61, 64, 
NOs Vey 74,075," (Oy 714 . 04. 0, 


Be Namelist AEROSO 


This namelist contains the assumed CCN profile (variable AERP, NZ 
values). Since in the present application the cloud module is 
turned off (IFCLOU=0) the values of AERP are irrelevant. 
Exemplary values are given below: 


AERP = 300,300,300,300,150,150,30,30,30,30,30,30,30,20,20,10 
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APPENDIX C 


Modifications of Gesima undertaken under Phase II 
under this project 
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APPENDIX C 


Modifications of Gesima undertaken under Phase II of this 
project. 


C.1 The MAIN program 


The new subroutine NUDEV (see Section C.4) is executed just 
before PRINTT at INCOUT intervals. 


Calls to subroutines WRTAP and RLUSE have been replaced by 
calls to new subroutines WRTAPF (see Section ) and RLUSEX (see 
Section C.8) 


C.2 Subroutine INITS 


The purpose of this subroutine is to initialize the 
variables of the surface layer module. In particular, the 
variable XSOIL contains for each grid cell the 8 surface 
parameters, computed as averages of values prescribed to several 
land use classes, weighted by the relative coverage of the given 
grid cell by those land types (see Section Al.2 of the last year 
report). | 


In place of 11 land use classes used in Phase I, the 8 
classes of ADOM (water, deciduous forest, coniferous forest, 
swamp, cultivated, grass, urban and desert) were introduced to 
Gesima. The surface parameters were assigned to these classes by 
directly assigning or, in some cases, averaging the values from 
the appropriate classes of the 19 used in the original Gesima. 


Reduction of the number of classes resulted in changes in 
dimensions of some variables (XWTYP, XSOIL etc.) and in the 
averaging formula. 


Analysis of the code of subroutine SURFW showed, that in 
some cases two XSOIL surface parameters were used in one 
statement. For example the product of the diffusion coefficient 
and the heat capacity XSOIL(i,j,1)*XSOIL(i,j,2) appeared in some 
statements. Thus, the product of two averages was used where an 
average of a product would be more appropriate. In order to 
correct this situation a new variable DIFCAP(NX,NY) - an average 
product of the diffusion coefficient and the heat capacity, 
weighted by the coverage of the grid cell by various land use 
classes is computed in INITS. 
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Similarly, the new variable AWFCX(NX,NY) is the averace 
product of field capacity and capillarity and DELTVX is the 
average reciprocal of field capacity multiplied by the time step 
DELTAT. 


A modification has been also introduced to the 
initialization of the wetness parameter SWET. For most grid cells 
it assumes, as before, the value obtained in the run of the 
one-dimensional Gesima preprocessor. For the cells containing 
more then 80% of water surface SWET is changed be equal to that 
water coverage (a number between 0.8 and 1.0). 


C.3 Subroutine INTNU 


Since the uniform cloudiness has been replaced in the 
subroutine RADSFC (see Section C.7) by a two-dimensional variable 
SKYCO, the interpolation in space of its values from the large 
scale low resolution data has been included in INTNU. 


Similarly, the interpolation in space of the mixing height 
XMIX, needed in the new version of MUGRAD (see Section C.5) has 
also been included. After interpolation from the large scale 
data, the Gesima terrain following height transformation is 
applied to XMIX. 


The interpolation in time of SKYCO and XMIX is performed in 
subroutine TDBCNU. 


C.4 Subroutine NUDEV 


NUDEV is a new auxiliary subroutine, computing for output 
purposes the deviations of the mesoscale fields of wind, 
temperature and humidity computed by Gesima from their smooth, 
low resolution values interpolated from the large scale data and 
used in the "nudging" scheme. 


The subroutine is called from the main program at INCOUT 
intervals. The three-dimensional deviation fields can be then 
written to output files in subroutine PRINTT. The vertical 
profiles of the minimum, maximum and mean deviations are written 
to the general output file (fort.06) in NUDEV. 
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C.5 Subroutine MUGRAD 


MUGRAD is the main subroutine of the boundary layer module 
of Gesima. It has been modified, to include the new convective 
PBL parameterization, described in Section 4 of this report. 


C.6 Subroutine PRINTT 


This subroutine writes three-dimensional Gesima variables to 
the output files. For details see Section A1.9 of the last year 
report. The subroutine required extensions to accommodate the 
deviation arrays computed in subroutine NUDEV (see Section C.4). 
Instead of adding 4 more output files and increasing the size of 
the program, the deviation variables replaced in PRINTT some of 
the less needed output fields. 


The temperature deviation is written out in place of the 
potential temperature to the unit IFOUTP(1). The deviations of u 
and v components replace the longwave radiation cooling and mass 
flux divergence in units IFOUTP(5) and IFOUTP(6), respectively. 
The humidity deviation from the large scale field replace the 
temperature deviation from the reference profile in IFOUTP(9). 


C.7 Subroutine RADSFC 


RADSFC computes the shortwave and longwave radiation fluxes 
at the ground surface. In the original Gesima two alternative 
approaches are used for parameterization of cloud effects on both 
solar and atmospheric radiation. For runs with the cloud module 
active, the "liquid water path", computed from the cloud liquid 
water content is used in RADSFC. If the cloud module is 
deactivated by setting the switch IFCLOU to 0 and the liquid 
water content is not available, Gesima uses the input cloud 
coverage parameter COV in a simplified set of radiation formulae. 


Gesima is able to simulate only the stratiform clouds for 
the horizontal resolution of 20km, so the liquid water path 
approach was deactivated in Phase I of this project, and COV used 
instead, even in runs with cloud module active. 


Extending the domain of Gesima made the use of a single 
input cloudiness parameter no longer reasonable. COV has been 
replaced in RADSFC by a two-dimensional array SKYCO, whose values 
are interpolated in time and space, from the large scale input 
data in subroutines INTNU and TDBCNU similarly to the 
meteorological fields used for "nudging" purposes. 
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C.8 Subroutine RLUSEX 


Since the structure of the input land use data has changed, 
the subroutine RLUSE, reading and interpolating the land use 
information under the Phase I has been replaced by a new 
subroutine RLUSEX. 


The interpolation of land use data is no longer necessary. 
RLUSEX reads the relative coverage of each grid cell by 8 classes 
of land use described in Section C.2 from eight input files: 

FOEL si cOrt<aL2, LOrt.235,.f0Gt.24, fort.35,. ,0Ort.36, fort.37 and 
fort-38: 


C.9 Subroutine RESTAR 


Statements reading the new variables AWFCX, DELTVX and 
DIFCAP, computed in the subroutine INITS (see Section C.2) and 
written to the restart file in SAFETY, have been added to RESTAR. 


C.10 Subroutine SAFETY 


Statements writing to the restart file the new variables 
AWFCX, DELTVX and DIFCAP, computed in the subroutine INITS (see 
Section C.2), have been added to SAFETY. 


C.11 Subroutine SURFW 


The new variables DIFCAP, AWFCX and DELTVX, computed in 
subroutine INITS (see Section C.2) replaced reciprocals or 
products of XSOILS variables, where appropriate. 


C.12 Subroutine TDBCNU 


Statements interpolating in time the new variables XMIX and 
SKYCO have been added to TDBCNU. 


C.13 Subroutine WRTAP 


The purpose of this subroutine is to write the 
three-dimensional fields generated by Gesima to an unformatted 
output file, for later use by ADOM. The version used under Phase 
I wrote out more data than were actually needed, in order to 
enable the changes in the Gesima/ADOM interface without the need 
of rerunning Gesima. After gaining more experience, some of those 


A LS 


C5 


redundant fields have been removed in Phase II in order to save 
disc space. 


C.14 Subroutine WRTAPF 


After transferring the system from the Cray computer of the 
Ontario Centre of Large Scale Computations to the Centre 
Informatique de Dorval (CID), Gesima is run on the SX-3 
supercomputer, while the Gesima/ADOM interface and ADOM use the 
CID Cray temporarily. This setup made the use of unformatted 
output file generated by WRTAP impossible. 


Until ADOM is moved to SX-3, subroutine WRTAP must be 
replaced by its new, temporary version WRTAPF, using formatted 
write statements, instead of unformatted ones. 


Since in the latest runs of Gesima the cloud module was 
disabled, the writing out of the cloud data have been removed 
from WRTAPF (and their use from the Gesima/ADOM interface). If 
the cloud module is reactivated, the removed statements in WRTAPF 
should be reactivated as well. 
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Gesima subroutines 


D1 
APPENDIX D 
Gesima Subroutines 


This Appendix provides a list of all Gesima subroutines with 
short descriptions of their purposes and the meaning of formal 
parameters. For more information on subroutines developed or 
substantially modified under this project, see the Phase I report 
(Niewiadomski and Shenfeld, 1991). The modifications made in 
Phase II are described in detail in Appendix C of this report. 


SUBROUTINE ADV(DC,C,ITYP) 


Purpose: Compute advective terms of a central variable C 
and add them to the increment DC 


Called from: ADVECN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 


Formal parameters: 


DC(0:NXP1,0:NYP1,0:NZP1) Increment array 
C(0:NXP1,0:NYP1,0:NZP1) Advected variable 
ETYP Switch depending of the type of 


adveced variable 


SUBROUTINE ADVECN 


Purpose: Call ADV for mass fluxes and temperature, 
perform auxiliary computations 


Called from: MAIN 
Subroutines called: ADV 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE BMUL(WL,WLS) 


Purpose: Multiply B by WLS to obtain WL. 
Auxiliary subroutine of the pressure solver. 
Called from: IGCG 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 


Formal parameters: WL(NVOL), WLS(NVOL) 
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SUBROUTINE BORLAN(C) 


Purpose: Apply Orlanski-type boundary conditions for 
variable C using phase velocity of temperature 


Comments : This subroutine was used in the orginal model 
for temperature, humidity and hydrometeor mixing 
ratios. After introduction of subroutine BORLAQ 
is used only for temperature. 


Called from: MAIN (if IFTINF.LE.0.OR.IFTOUT.LE.0) 
Subroutines called: None 

Origin: GKSS 

Main modifications: None 

Formal parameters: C(0:NXP1,0:NYP1,0:NZP1) 


SUBROUTINE BORLAQ(C) 


Purpose: Apply Orlanski-type boundary conditions for 
variable C using phase velocity of humidity. 


Comments: This is a modified subroutine BORLAN with phase 
velocities VPHT* replaced by VPHQ*. Applied for 
humidity and hydrometeor mixing ratios. 


Called from: CLOUDP, HUMID (if IFQINF.LE.0.OR.IFQOUT.LE.0) 


Subroutines called: None 
Origin: MEP 
Formal parameters: C(O0:NXP1,0:NYP1,0:NZP1) 


SUBROUTINE BOTTOM 


Purpose: Compute u*, T* and q* if surface module is not 
called 
Comments: Not used, since the surface module is always 


active in these simulations. 


Called from: MAIN (if IFUBOT.EQ.3.AND.IFSURF.EQ.0) 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 


Formal parameters: None 
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SUBROUTINE CBOUND(C,IFI) 
Purpose: Apply non-Orlanski lateral boundary conditions 


to variable C (temperature, humidity or 
hydrometeor mixing ratio) 


Comments : IFI is equal to IFTINF for temperature or 
to IFQINF for other variables 
Called from: MAIN, CLOUDP, HUMID 
Subroutines called: None 
Origin: GKSS 
Main modifications: switch IFTINF is replaced by IFI which 


has been added as a formal parameter. 


Formal parameters: 

C(0:NXP1,0:NYP1,0:NZP1) the variable processed 

IFI switch definig the type of boundary 
conditions 


SUBROUTINE CENTUP 


Purpose: Update the ‘central variables’ (wind components 
and temperature) by adding their increments 


Comments : Humidity and hydrometeors are updated in HUMID and 
CLOUDP, respectively 

Called from: MAIN 

Subroutines called: None 

Origin: GKSS 

Main modifications: Nudging terms added 

Formal parameters: None 


SUBROUTINE CLOUDP 


Purpose: Compute transformations and transport of 
hydrometeors. 
Comments : Driving subroutine of the cloud module. 


Not active in recent applications. 
Called from: MAIN (if IFCLOU.EQ.1) 
Subroutines called: TRANSP, DIFFC, TEMP, BORLAQ, CBOUND, 


QLIRTB 
Origin: GKSS 
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Main modifications: Calls to BORLAN replaced by ones to 
BORLAQ 
IFQINF included in CALL CBOUND 
statements. 
Formal parameters: None 


SUBROUTINE CONCMC(XCMC,YCMC,ALAT,ALON) 


Purpose: Convert CMC coordinates to lat-lon 

Comments : Based on the program CONVCMC received from OME. 
Used in older versions (which use RLUSE) only 

Called from: RLUSE 

Origin: MEP 

Formal parameters: 

XCMC, YCMC CMC coordinates (input) 

ALAT, ALON latitude, longitude (output) 


SUBROUTINE CONST 


Purpose: Initialize various physical constants 
Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE CORIOL 


Purpose: Compute the Coriolis force terms 

Called from: MAIN (if IFCORI.GT.0) 

Subroutines called: None 

Origin: GKSS 

Main modifications: Uniform geostrophic wind replaced by a 


3-D array 


Formal parameters: None 
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SUBROUTINE CUMU 
Purpose: Compute the tendency terms for hydrometeors 
Comments: This is the main subroutine of the cloud module. 


Not active in recent applications. 


Called from: CLOUDP 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE DEFRAD 


Purpose: Define some constants for the radiation module 
Called from: MAIN (if IFRAD.EQ.1) at the beginning of each run 


Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE DIFFC(DC,C,M1T,M3T,ITYP) 


Purpose: compute the diffusive source term for a centered 
variable C 
Called from: CLOUDP, DIFFIL, HUMID (if IFDIFF.NE.0) 


Subroutines called: IMPDI2 
Origin: GKSS 
Main modifications: None 


Formal parameters: 


C(0:NXP1,0:NYP1,0:NZP1) variable processed 
DC(0:NXP1,0:NYP1,0:NZP1) increment array 
MIT(0O:NXP1,0:NYP1,0:NZP1) horizontal diffusivity 
M3T(0:NXP1,0:NYP1,0:NZP1) vertical diffusivity 
TYP passed to IMPDI2 


SUBROUTINE DIFFIL 


Purpose: Call DIFFC and FILTER for wind and temperature. 
Convert temperature to real value before that and 
back to deviation form after. 

Called from: MAIN (if IFDIFF.NE.0.OR.IFFILT.NE.0) 
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Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE DISPLU 
Purpose: Compute velocities from fluxes 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE FBOUND 
Purpose: Compute non-Orlanski boundary terms for wind 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE FGRAD(C) 
Purpose: Compute pressure gradient and update fluxes 


to conserve mass. 


Called from: 
Subroutines called: 
Origin: 

Main modifications: 


Formal parameters: 
C(0:NXP1,0:NYP1,0:NZP1) 


SUBROUTINE FILTER(DC1D,C) 


MAIN 
None 
GKSS 
None 


variable processed. Actually XPRES 
only. 


Purpose: Apply 3 or 5-point filter to variable C 
Called from: DIFFIL, HUMID (if IFFILT.NE.0) 


Subroutines called: 
Origin: 


None 
GKSS 
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Main modifications: None 
Formal parameters: None 
SUBROUTINE FLUCOR 
Purpose: Compute new inflow rate, check balance and adjust 

outflow. 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE FLUXFI 
Purpose: Find final fluxes by taking the mean of old and 
corrector values 

Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE FLUXUP 
Purpose: Compute fluxes from velocities 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE FORLAN 
Purpose: Apply Orlanski boundary conditions to fluxes 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 


Formal parameters: None 
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SUBROUTINE HUMID 

Purpose: Solve the humidity equation. 

Called from: MAIN (if IFHUMI.EQ.1) 

Subroutines called: TRANSP, DIFFC, FILTER, BORLAQ, 
CBOUND, QVTOB 

Origin: GKSS 

Main modifications: Call to BORLAN replaced by one to 
BORLAQ. 
IFQINF added as formal parameter of 
CBOUND. 


Nudging terms added 


Formal parameters: None 


SUBROUTINE HYHORI 


Purpose: Compute horizontal hydrostatic pressure gradient 
Called from: MAIN (if IFHYGR.EQ.1) 


Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE HYVERT 


Purpose: Compute vertical hydrostatic pressure gradient 
Called from: MAIN (if IFHYGR.EQ.1 and IFSEMI.NE.0) 


Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE IBOUND 
Purpose: Initialize some surface parameters if IFSURF = 0 


Comments: Not used in this application, since the surface 
module is always active (IFSURF = 1) 


Called from: MAIN (if IFSURF.EQ.0 and IFHOTS.EQ.0) 


Subroutines called: None 
Origin: GKSS 
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Main modifications: None 
Formal parameters: None 
SUBROUTINE IGCG 
Purpose: Solves a linear system of equations by 

preconditioned conjugate gradient method 

Comments : Main subroutine of the pressure solver 
Called from: MAIN 
Subroutines called: BMUL, LIMUL, UIMUL 
Origin: GKSS 
Main modifications: Some informative output removed 
Formal parameters: None 


SUBROUTINE IJKPOS(INDEX,I,J,K) 


Purpose: Calculate indices (i,j,k) from 1-D index of a 
centered variable 

Comments: Auxiliary subroutine, not used in this version 

Called from: not used in this version 

Subroutines called: None 

Origin: GKSS 

Main modifications: None 


Formal parameterers: 
INDEX index of a 1-D array (input) 
ŒUrK indices of the equivalent 3-D array (output) 


SUBROUTINE IMBOUY 


Purpose: Computes buoyancy terms 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE IMPDI2(DC1D,C1D,IM3T2D,ITYP) 


Purpose: Compute implicit vertical diffusion 
Called from: DIFFC 
Subroutines called: LUSOL2 


Origin: GKSS 
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Main modifications: Not used part for ITYP=4 removed 
Formal parameters: 
C1D(NVOL) variable processed 
DC1D(NVOL) incerment array 
IM3T2D(NH2,0:NZP1) diffusivity array 
MNS Type of variable processed: 

ITYP = 1 for temperature, moisture, hydrometeors 

ITYP = 2 for horizontal velocities 

ITYP = 3 for vertical velocity 

ITYP = 4 for passive tracer 


SUBROUTINE INIFLA 


Purpose: Initialize fluxes according to initial velocity 
fields 
Comments: based on subroutine FLUXUP 


replaces subroutine INITFL 
Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: None 
Origin: MEP 
Formal parameters: None 


SUBROUTINE INITCL 


Purpose: Initialize hydrometeor fields and some contstants 
of the cloud module 


Comments: Hydrometeors set to 0 for initial runs only 


Called from: MAIN at the beginning of each run 


Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE INITFL 
Purpose: Initialize inflow flux profiles 


Comments: Not used in this version, replaced by INIFLA 
Called from: Not used in this version, replaced by INIFLA 


Subroutines called: None 
OrExgini: GKSS 
Main modifications: None 


Formal parameters: None 
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SUBROUTINE INITFT 


Purpose: Adjust boundaries for flux-conservation and 
compute total flux 
Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE INITNU 


Purpose: Initialize interpolation coefficients for nudging 
and boundary conditions. 
Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: None 
Origin: MEP 
Formal parameters: None 


SUBROUTINE INITS 


Purpose: Initialize surface model constants and variables 

Called from: MAIN (if IFSURF.EQ.1) at the beginning of each 
run 

Subroutines called: None 

Origin: GKSS 

Main modifications: SOILS redefined, XSOIL introduced. 


DIFCAP, AWFCX, DELTVX variables 
introduced in Phase II, see 
Appendix C) 


Formal parameters: None 


SUBROUTINE INITUR 


Purpose: Initialize constants of the turbulence module 
Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: KBOUND 
Origin: GKSS 
Main modifications: None 


Formal parameters: None 
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SUBROUTINE INTBC 
Purpose: Interpolate time dependent boundary conditions 
Comments: Mostly obsolete, since for most variables t.d.b.c. 


have been replaced by nudging. 


Called from: MAIN at the beginning of each run and every full 


hour 
Subroutines called: None 
Origin: MEP 
Formal parameters: None 


SUBROUTINE INTINI 


Purpose: Initialize interpolation coefficients for 
boundary conditions 


Comments: Used in pre-August 1991 versions only 


Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: None 
Origin: MEP 
Formal parameters: None 


SUBROUTINE INTNU 


Purpose: Interpolate nudging fields, sky coverage and 
mixed layer height. 


Comments: Modified in Phase II, see Appendix C. 


Called from: RDBCI at the beginning of an initial run and MAIN 
at the beginning of a restart run and every full 


hour 
Subroutines called: None 
Origin: MEP 
Formal parameters: None 


SUBROUTINE INZRTP 


Purpose: Initialize reference density, Orlanski phase 
velocities and some other variables. 
Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


« à £Æ \ 


D13 
Subroutines called: None 
Origin: GKSS 
Main modifications: Phase velocities for humidity type 


variables (VPHQE etc.) added. 


Formal parameters: None 


SUBROUTINE KBOUND 


Purpose: Apply boundary conditions for turbulent 
diffusion coefficients 

Comments: Active only if IFDIFF.GE.2 

Called from: INITUR, MUGRAD 

Subroutines called: None 

Origin: GKSS 

Main modifications: None 

Formal parameters: None 


SUBROUTINE LIMUL(WLS,WL) 


Purpose: Multiply the inverse of L with WL to obtain WLS 
Called from: IGCG 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 


Formal parameters: 
WL (NVOL) operand array 
WLS ( NVOL) product array 


SUBROUTINE LUSOL2(A,B,C,F,X,N,NH2) 


Purpose: solve a tridiagonal matrix using Thomas algorithm 
Called from: IMPDI2, SURFW 

Subroutines called: None 

Origin: GKSS 

Main modifications: None 


Formal parameters: 

A(NH2,0:N),B(NH2,0:N),C(NH2,0:N),F(NH2,0:N) coefficient arrays 
X(NH2,0:N) processed array 
N,NH2 array dimensions 
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SUBROUTINE LWRAD 
Purpose: Compute longwave radiation fluxes and flux 
divergence 
Comments : Not used in this application (IFRAD set to 0) 
Called from: MAIN (if IFRAD.EQ.1) 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE MOMFLU 
Purpose: Compute momentum flux for output purposes 
Called from: MAIN at INCOUT intervals 
Subroutines called: None 
Origin: GKSS 
Main modifications: 3-D geostrophic wind introduced 
Formal parameters: None 
SUBROUTINE MOMFNO 
Purpose: Compute analytic hydrostatic momentum flux 


for scaling 


Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE MONITO 


Purpose: Compute and print max and min values of some 
parameters 

Comments: Used to monitor propper behaviour of the model 

Called from: MAIN at IFMONI intervals 

Subroutines called: None 

Origin: GKSS 

Main modifications: None 


Formal parameters: None 


SUBROUTINE MUGRAD 
Purpose: 
Comments: 
Called from: 
Subroutines called: 
Origin: 


Main modifications: 
Formal parameters: 


SUBROUTINE NEWOLD 


Purpose: 


Called from: 


Subroutines called: 
Origin: 


Main modifications: 
Formal parameters: 


SUBROUTINE NUDEV 


Purpose: 
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Compute eddy diffusivity coefficients 


New convective PBL parameterization introduced 


MAIN (if IFDIFF.EQ.2) 
STRATI, KBOUND 

GKSS 

See Appendix C 

None 


Replace ‘old’ velocities and temperature with 
the ‘new’ ones. 


MAIN at the end of each time step 


None 
GKSS 


None 
None 


Compute deviations from the nudging fields for 


output purposes. 
Compute and print max, min and average values 
of velocity, temperature and humidity 


Comments : 


Called from: 
Subroutines called: 
Origin: 

Formal parameters: 


SUBROUTINE NUDGING 


Purpose: Compute the 
Called from: 

Subroutines called: 
Origin: 

Formal parameters: 


Developed in Phase II, see Appendix C. 


MAIN at INCOUT intervals 
None 

MEP 

None 


nudging terms 


MAIN 
None 
MEP 

None 
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SUBROUTINE ORLAN 


Purpose: Set phase velocities equal to outflow velocities 
and determine the Courant number 


Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: Separate (from temperature) phase 


velocities for humidity and 
hydrometeors introduced 
Formal parameters: None 


SUBROUTINE OUTERF 


Purpose: Adjust outer values for all fluxes 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE PHYDRO(C1D) 


Purpose: Compute the hydrostatic pressure field 
Called from: MAIN (if IFHYUP.EQ.1) 

Subroutines called: None 

Origin: GKSS 

Main modifications: None 

Formal parameters: 

C1D(NVOL) Temperature deviation array 


SUBROUTINE PLHS 


Purpose: Compute the left hand side of the Poisson equation 
for pressure 
Comments: Part of the pressure solver package 


Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: None 
Origin: GKSS 
Main modifications: None 


Formal parameters: None 
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SUBROUTINE PRECG 
Purpose: Compute the right hand side of the Poisson 
equation for pressure 
Comments: Part of the pressure solver package 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE PRIGEO 
Purpose: Print geometry of grid 
Comments: Auxiliary subroutine, not used in this version 
Called from: not used in this version 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE PRINTT 
Purpose: Write 3-D arrays of specified parameters to output 
files 
Comments: Main output subroutine of the model 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: Some formats changed. 
Some fieds redefined in Phase II, 
see Appendix C. 
Formal parameters: None 


SUBROUTINE QLIRTB 


Purpose: Assign top and bottom boundary conditions for 
hydrometeors 

Called from: CLOUDP 

Subroutines called: None 

Origin: GKSS 

Main modifications: None 


Formal parameters: None 
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SUBROUTINE QVTOB 
Purpose: Assign top and bottom boundary conditions for 
humidity 
Called from: HUMID 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE RADSFC 
Purpose: Compute surface radiation fluxes 
Comments : For IFRAD=0 (as is the case for this application) 


this subroutine parameterizes the longwave 
radiation, which for IFRAD=1 is computed in LWRAD 


Modified in Phase II, see Appendix C. 


Called from: MAIN (if IFSURF.EQ.1) 
Subroutines called: None 

Origin: GKSS 

Main modifications: See Appendix C. 
Formal parameters: None 


SUBROUTINE RDBCI 


Purpose: Read boundary conditions and nudging data. 
Initialize some variables. 


Comments: A version of RDBC for the initial runs. 
Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: INTNU 
Origin: MEP 
Formal parameters: None 


SUBROUTINE RDBC 
Purpose: Read boundary conditions and nudging data. 


Called from: MAIN at the beginning of a restart run and every 
full hour 


Subroutines called: None 
Origin: MEP 
Formal parameters: None 
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SUBROUTINE RDTAPE 
Purpose: Read initial profiles from the output of 1-D run 
Comments: — Most parameters overwritten in RDBCI 


Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: None 

Origin: GKSS 

Main modifications: Structure of 1-D output and the 
read statements slightly changed. 

Formal parameters: None 


SUBROUTINE READD 


Purpose: Read the input file and initialize some constants 
Comments: For the detailed structure of the input file see 
Appendix B 


Called from: MAIN at the beginning of each run 


Subroutines called: None 
Origin: GKSS 
Main modifications: Some namelists redefined. 


Reading topography data moved to 
RDPROF (for IFTOPO = 0 only). 
Grid data read in from input 
redefined inside the subroutine. 


Formal parameters: None 


SUBROUTINE RLUSE 


Purpose: Read and interpolate land use data. 

Comments : Used in older versions of the model. If the 
interpolated land use data are available, RLUSEX 
is called instead. 


Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 
Subroutines called: CONCMC 


Origin: MEP 
Formal parameters: None 
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SUBROUTINE RLUSEX 
Purpose: Read land use data. 
Comments: Developed in Phase II to replace RLUSE if the 


interpolated land use data are available. Each 
category is read from a different file.The 
subroutine is domain-dependent. For details see 
Appendix C. 


Called from: MAIN (if IFHOTS.EQ.0 - initial run only) 


Subroutines called: None 
Origin: MEP 
Formal parameters: None 


SUBROUTINE RESTAR 
Purpose: Read the unformatted restart file 


Called from: MAIN (if IFHOTS.EQ.1 - at the beginning of a 
restart run) 


Subroutines called: None 
Origin: GKSS 
Main modifications: See Appendix C 
Formal parameters: None 


SUBROUTINE RTOPO 
Purpose: Read topography data 


Called from: MAIN (if IFTOPO.EQ.0 and IFHOTS.EQ.0) at the 
beginning of an initial run only. 


Subroutines called: None 
Origin: MEP 
Formal parameters: None 


SUBROUTINE SAFETY (NTAPE) 

Purpose: Write an unformatted restart file 

Comments: A single NTAPE value and single calls to SAFETY 
used in this version, instead of the double calls 


in original Gesima. 


Called from: MAIN at INCHOT intervals 
Subroutines called: None 
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Origin: GKSS 
Main modifications: See Appendix C 
Formal parameters: None 
SUBROUTINE SETIND 
Purpose: Set indicator arrays for inflow/outflow 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE SHIFTT 
Purpose: Define shift coefficients 


Compute velocities from fluxes at the beginning 
of a corrector step. 


Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE SPRINT 


Purpose: Write some surface module parameters to an output 
file 


Called from: MAIN at incout intervals (if IFSURF.EQ.1) 


Subroutines called: None 
Origin: GKSS 
Main modifications: Some parts removed 
Formal parameters: None 


SUBROUTINE STAU 


Purpose: Compute auxiliary variable TAU1D 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 


Main modifications: None 
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Formal parameters: None 
SUBROUTINE STRATI 
Purpose: Compute auxiliary terms for diffusivity 
parameterization 
Called from: MUGRAD 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE SURFW 
Purpose: Compute the surface energy balance 
Comments: Main subroutine of the surface module 
Called from: MAIN (if IFSURF.EQ.1) 
Subroutines called: LUSOL2 
Origin: GKSS 
Main modifications: SOILS replaced by XSOIL. 


DIFCAP, AWFCX and DELTVX 
introduced, see Appendix C. 


Formal parameters: None 


SUBROUTINE TDBCNU 


Purpose: Assume new values of nudging variables, cloudiness 
and mixed layer height. 
Assume new inflow boundary conditions. 


Comments : The scheme defining the boundary conditions for 
hydrometeors not tested sufficiently (but the 
cloud module not active in recent simulations). 
Interpolation of cloudiness and mixing height 
introduced in Phase II. 


Called from: MAIN 
Subroutines called: None 
Origin: MEP 


Formal parameters: None 
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SUBROUTINE TEMP 


Purpose: Convert temperature deviation TNC to actual 
absolute temperature. 


Comments : Used mainly for output purposes 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE TIMSTP 


Purpose: Compute DELTAT from numerical stability 
criteria for runs with variable time step 


Comments : Not tested and not used in this application, 
(IFTIM set to 0 - constant time steps) 


Called from: MAIN (if IFTIM.EQ.1) 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 


SUBROUTINE TRANSP(DC1D,C1D,W) 


Purpose: Compute advective transport of passive components 
using Smolarkiewicz scheme 


Called from: CLOUDP, HUMID 
Subroutines called: None 

Origin: GKSS 

Main modifications: None 

Formal parameters: 

C1D(NVOL) processed variable 
DC1D(NVOL) increment array 
W(0:NXP1,0:NYP1,0:NZP1) sedimentation velocity 


SUBROUTINE TTOP 


Purpose: Top boundary conditions for temperature 
Called from: MAIN 
Subroutines called: None 


Origin: GKSS 
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Main modifications: None 
Formal parameters: None 
SUBROUTINE UIMUL(WLS,WL) 
Purpose: Multiply the inverse of U by WL to obtain WLS 
Comments : Part of the pressure solver package 
Called from: IGCG 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: 
WL (NVOL) operand array 
WLS (NVOL) product array 
SUBROUTINE UVTBOT 
Purpose: Lower boundary conditions for U, V , W, T 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
SUBROUTINE VERIFY 
Purpose: Print some input/restart parameters 
Comments: Auxiliary subroutine 


Called from: MAIN at the beginning of each run 


Subroutines called: None 
Origin: GKSS 
Main modifications: Some output removed 
Formal parameters: None 


SUBROUTINE WRTAP(NTAPE) 
Purpose: Write binary output data for ADOM purposes 
Comments: Based on subroutine SAFETY. 


Some fields removed in Phase II - see Appendix C. 
Temporarily replaced by WRTAPFX - see Appendix C. 


Called from: MAIN at INTAP intervals 
Subroutines called: None 
Origin: MEP 


Formal parameters: None 
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SUBROUTINE WRTAPF (NTAPE) 
Purpose: Write formatted output data for ADOM purposes 
Comments: Temporary replacement of subroutine WRTAP, see 
Appendix C. 
Called from: MAIN at INTAP intervals 
Subroutines called: None 
Origin: MEP 
Formal parameters: None 
SUBROUTINE WTOPBO 
Purpose: Update top and bottom b. c. for vertical flux 
Called from: MAIN 
Subroutines called: None 
Origin: GKSS 
Main modifications: None 
Formal parameters: None 
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